function DC_spike_xcorr(varargin)
% cross correlation of spike trains -
% version to work with dynamic clamp data from pairs of cells.
% 5/14/2006 P. Manis...
% assumes data are in channels 1 and 3 of each selected dataset, and
% computes (for each trace) an xcorr, and sums the xcorrs through the
% dataset and normalizes. 
% if no spikes have been analyzed, then you'll need to to that first.
% Results are stored in XCRESULT... not in CONTROL. At present, this must
% be manually saved.
% If multiple database entries are selected, the routine will compute the
% ccross correlations between each unique pair of spike trains and store
% them in XCRESULT.
% Note: an estimate of the cross-correlation significance is obtained by
% shuffling the intervals in one of the spike trains, while retaining the
% first spike latency and the interspike interval distribution (e.g.,
% interval order is shuffled). The number of iterations is controlled by
% nshuffle in the template_go routine. large nshuffle can take a VERY long
% time to compute however.
%%
% TESTMODES:
% I provide three test modes for this routine. 
% TESTMODE = 0 - runs on the data selection
% TESTMODE = 1 - Two sets of uncorrelated spike trains are generated.
% based on:
% basic template file for data analysis routines.
% gets the current datac window selection, and does the analysis
% for all selected record groups/blocks
% automatically refreshes the data arrays from the disk file if necessary
% and also shows how to generate a figure.
% 8/23/05 P. Manis

% the main routine is just a wrapper to handle the multiple selection and
% catch execution errors.
%
% Note: when testing, set the testflag to 1 so that we DO NOT do try/catch.
% you need to do this to catch the errors in the _go routine.
%
global CONTROL
global XCRESULT

testflag = 1;

if(nargin >= 1)
    testmode = varargin{1};
else
testmode = 0; % see above for testmode descriptions.
end;
if(testmode == 0)
    sf = getmainselection; % read the list of selected traces
else
    sf = [1,2]; % fake it.
end;


if(sf > 0)
    pflag = getplotflag; % read the plot flag (used in the routine template2)
    XCRESULT = []; % initialize with every run to minimize confusion.
    QueMessage('Spike train xcorr analysis routine', 1); % clear the que
    k = 1;
    for i = 1:length(sf)
        if(testflag == 0)
            try
                template_go(sf(i), k, testmode, pflag); % between selected pairs of datasets...
            catch
                if(testmode == 0)
                    fprintf(1, 'Spike train xcorr Analysis Failed execution for sf = %d\n  ***[File: %s, records: %d - %d ]\n', ...
                        sf(i), CONTROL(sf(i)).filename, CONTROL(sf(i)).recbeg, CONTROL(sf(i)).recend);
                else
                    fprintf(1, 'Spike train xcorr TESTMODE %d failed\n', testmode);
                end;
            end;
        else
            template_go(sf(i), k, testmode, 0); % also turn off plot flag when testing
        end;

    end;
    k = k + 1;

end;



%--------------The real analysis routine------------------

function template_go(sf1, sf2, k, testmode, pflag)

% access to all of the data and parameters
global  DFILE CONTROL ALLCH XCRESULT
%
% paramteres:
minlat = 100;
maxlat = 900;
tmaxlag = 500;
h = findobj('tag', 'simcorrfig');
if(isempty(h))
    h = figure;
    set(h, 'tag', 'simcorrfig');
else
    figure(h);
    clf;
end;
hax1 = subplot('Position', [0.1 0.67 0.8 0.25]);
hax1a = subplot('Position', [0.1,0.37, 0.8, 0.25]);
hax2 = subplot('Position', [0.1 0.22 0.8 0.10]);
hax3 = subplot('Position', [0.1 0.07 0.8 0.10]);
htitle = subplot('Position', [0.1 0.93 0.8 0.07]);
axis([0,1,0,1])
axis('off')
if(testmode == 0)
fn1 = CONTROL(sf1).filename;
recb1 = CONTROL(sf1).recbeg;
rece1 = CONTROL(sf1).recend;
comm1 = CONTROL(sf1).Commentary;
fn2 = CONTROL(sf2).filename;
recb2 = CONTROL(sf2).recbeg;
rece2 = CONTROL(sf2).recend;
comm2 = CONTROL(sf2).Commentary;

text(0,0.8126,sprintf('%-12s R[%d:%d] (%s) VS %-12s R[%d:%d] (%s)', ...
    CONTROL(sf1).filename, CONTROL(sf1).recbeg, CONTROL(sf1).recend, CONTROL(sf1).Commentary,...
    CONTROL(sf2).filename, CONTROL(sf2).recbeg, CONTROL(sf2).recend, CONTROL(sf1).Commentary...
    ), 'Fontsize', 8);


% Only compute if the spike trains have already been computed.

if(isempty(CONTROL(sf1).spike) | isempty(CONTROL(sf2).spike))
    QueMessage(sprintf('Missing spike analysis in pair %d, %d', sf1, sf2), 1);
    return;
end;

a = CONTROL(sf1).spike;
b = CONTROL(sf2).spike;
else
text(0,0.8126,sprintf('TESTMODE %d',testmode), 'Fontsize', 8);
fn1 = 'testcell1';
recb1 = 1;
rece1 = 50;
comm1 = '';
fn2 = 'testcell2';
recb2 = 1;
rece2 = 50;
comm2 = '';
[a,b] = makeruns(testmode, 0);


end;

% compute cross correlation of the spike trains here
% simxcorr2 does this as a sum of the cross-correlations between trains in
% the first two arguments, on a trial-by-trial basis (as if simultaneously
% recorded).
% Plot into the axis pointed by hax1.
[cx, lags] = simxcorr2(a.latency, b.latency, hax1, [minlat maxlat], tmaxlag);
yl = get(gca, 'YLim');
set(gca, 'YLim', [0 yl(2)]);
drawnow
%
% now compute a shuffled xcorr. Shuffle the 'a' spikes in each trial,
% maintaining the same isi distribution, and first spike latency, but all
% remaining spikes are repositioned.
%
nsa = 0; nsb = 0; tlag = []; mcx2 = []; scx2 = [];
nshuffle = 5; % number of different shuffles to average. Setting nshuffle to 0 skips this step
if(nshuffle > 0)
    for m = 1:nshuffle
        y={};
        for i = 1:length(a.latency)
            x = a.latency{i}; % get the list
            nsa = nsa + length(find(a.latency{i} > minlat & a.latency{i} < maxlat));
            nsb = nsb + length(find(b.latency{i} > minlat & b.latency{i} < maxlat));
            intvls = diff(x);
            j = randperm(length(intvls)); % randomly reorder the list
%            intvls = intvls(j); % reorder the intervals
            xy(1) = x(1); % keep the first spike latency
            xy(2:length(intvls)+1) = xy(1)+cumsum(intvls(j)); % turn intervals back into latencies
            %for k = 1:length(intvls)
            %    xy(k+1) = xy(k)+intvls(k); % and put into the y array
            %end;
            y{i} = xy;

        end;
        [cxs, lagshuffle, na, nb, dt] = simxcorr2(y, b.latency, hax2, [minlat maxlat], tmaxlag); % compute xcorr again...
        drawnow
        if(m == 1)
            cx2 = cxs;
        else
            cx2 = cx2 + cxs;
        end;
    end;
    %pause(2);
    %QueMessage('paused');
    cxshuffle = cx2/nshuffle;
    nsa = nsa/nshuffle;
    nsb = nsb/nshuffle;
    axes(hax2);
    bar(lagshuffle, cxshuffle);
    set(gca, 'Xlim', [-tmaxlag tmaxlag]);

    yl = get(gca, 'YLim');
    set(gca, 'YLim', [0 yl(2)]);
    set(gca, 'Xlim', [-100 100]);

    % plot the difference correlation (cross correlation minus shuffled average)
    axes(hax1a);
    plot(lagshuffle, cx-cxshuffle);
    set(gca, 'Xlim', [-100 100]);

    % compute some limits from the data in the shuffle.
    mcx2 = mean(cxshuffle);
    scx2 = std(cxshuffle);
    axes(hax1);
    hold on;
    tlag = get(gca, 'Xlim');
    plot(tlag, [mcx2 mcx2], 'r-');
    plot(tlag, [mcx2-scx2 mcx2-scx2], 'r--');
    plot(tlag, [mcx2+scx2 mcx2+scx2], 'r--');
    nsa = nsa / length(a.latency);
    nsb = nsb / length(b.latency); % spikes per trial.
    nsa = 1000*nsa / (maxlat - minlat); % convert to average rate over the window
    nsb = 1000*nsb / (maxlat - minlat);
    plot(tlag, [nsa nsa], 'g-');
    plot(tlag, [nsb nsb], 'c-');
    set(gca, 'Xlim', [-100 100]);
    drawnow
else
cxshuffle = [];
lagshuffle = [];
end;

XCRESULT(k).R1 = sf1;
XCRESULT(k).F1 = sprintf('%s[%d]', fn1, recb1);
XCRESULT(k).R2 = sf2;
XCRESULT(k).F2 = sprintf('%s[%d]', fn2, recb2);
XCRESULT(k).cx = cx; % cross correlation
XCRESULT(k).lags = lags; % lag times for correlation
XCRESULT(k).cxshuffle = cxshuffle; % shuffle cx for statistical comparison
XCRESULT(k).lagshuffle = lagshuffle; % lags for shuffle
XCRESULT(k).tlag = tlag; % time base for mean and stdev of shuffle values.
XCRESULT(k).meancx = mcx2; % mean and stdev of shuffle values.
XCRESULT(k).stdcx = scx2;

currentflag = 0; % set this to 1 to cross-correlate the current traces also -
% requires access to data and reading the data directly to obtain the
% current traces.


if(currentflag == 1)
    if(isempty(XCRESULT) | length(XCRESULT) < sf1 | ...
            isempty(XCRESULT(sf1)) | isempty(find(XCRESULT(sf1).R2 == sf2)) ) % conditions requiring calculation of cross correlation of current traces
        % set some control parameters
        isrc = 1;
        if(length(XCRESULT) < sf1)
            XCRESULT(sf1).R2=[];
            XCRESULT(sf1).xcorr=[];
            XCRESULT(sf1).xcorrlag=[];
        end;
        if(isfield(XCRESULT, 'R2'))
            nxc = length(XCRESULT(sf1).R2)+1; % indentify position to store result into
        else
            nxc = 1;
            XCRESULT(sf1).R2=[];
        end;
        j = find(XCRESULT(sf1).R2 == sf2);
        if(~isempty(j))
            nxc = j(1); % will overwrite the first incidence of the same match.
        end;
        % this next routine reads all of the data and information we need,
        % including calculating some standard stuff that is put into DPAR.
        % you will want DPAR.time especially - it is the time base that goes
        % with the data.

        [DFILE1, DPAR1, err] = ana_setup(DFILE, sf1);
        if(err ~= 0)
            return;
        end;

        if(~isempty(ALLCH))
            %    VOLTAGE = ALLCH{isrc};
            curr1 = ALLCH{isrc+1};
        end;

        [sfile1, df1] = block_info(sf1);

        [DFILE2, DPAR2, err] = ana_setup(DFILE, sf2);
        if(err ~= 0)
            return;
        end;
        [sfile2, df2] = block_info(sf2);

        if(~isempty(ALLCH))
            %    VOLTAGE = ALLCH{isrc};
            curr2 = ALLCH{isrc+1};
        end;
        % assumption: data in both traces sampled at the same rate...
        tmin = floor(minlat/DPAR1.rate(1));
        if(tmin <= 0)
            tmin = 1;
        end;

        tmax = floor(maxlat/DPAR1.rate(1));
        if(tmax > size(curr1,2))
            tmax = size(curr1, 2);
        end;
        maxlag = floor(tmaxlag/DPAR1.rate(1));
        [ixc, lags] = xcorr(curr1(1,tmin:tmax)-mean(curr1(1,tmin:tmax)),...
            curr2(1,tmin:tmax)-mean(curr1(1,tmin:tmax)), maxlag, 'coeff');
        XCRESULT(sf1).R2(nxc) = sf2;
        XCRESULT(sf1).xcorr{nxc} = ixc;
        XCRESULT(sf1).xcorrlag{nxc} = lags*DPAR2.rate(1);

        axes (hax3);
        nxc = find(XCRESULT(sf1).R2 == sf2);
        bar(XCRESULT(sf1).xcorrlag{nxc}, XCRESULT(sf1).xcorr{nxc});
        set(gca, 'XLim', [-tmaxlag +tmaxlag]);
        set(gca, 'Ylim', [0 1]);
        drawnow

    end;
end;
hfig = findobj('tag', 'simcorrfig');
if (getplotflag > 0)
    print (hfig, '-dps2');
end

QueMessage(sprintf('Spike xcorr analysis %d-%d done', sf1, sf2));



%----------------
% makeruns : generate two fake data sets for testing the routine above.

function [a, b] = makeruns(testmode, pflag)

start = 100;
baseper = 20; % 50 Hz.
stimdur = 1000;
ntestrep = 50;
a=[];
b=[];
switch(testmode)
    case 1 % exactly synchronized, and regular
        for i = 1:ntestrep
            a.latency{i} = [start:baseper:stimdur];
            b.latency{i} = [start:baseper:stimdur];
        end;

    case 2 % one is regular, the other has random jitter
        for i = 1:ntestrep
            X = [start:baseper:stimdur];
            a.latency{i} = X;
            b.latency{i} = X + 5*randn(size(X)); % normal random distrubiton on b dataset
        end;

    case 3 % both are regular, but both have independent random jitter
        for i = 1:ntestrep
            X = [start:baseper:stimdur];
            a.latency{i} = X + 5*randn(size(X));
            b.latency{i} = X + 5*randn(size(X)); % normal random distrubiton on b dataset
        end;

    case 4 % both are random and independent
        for i = 1:ntestrep
            rfac = 2;
            X=randpois(baseper*rfac, floor((stimdur-start)/(baseper*rfac))); % 1 adds fake refractory period
            a.latency{i} = start + cumsum(X); % corresponds to about 50 s/s
            X=randpois(baseper*rfac, floor((stimdur-start)/(baseper*rfac))); % separate process
            b.latency{i} = start + cumsum(X); % corresponds to about 50 s/s

        end;
        ddur = 100;

    otherwise
        return;
end;

if(pflag) % then try plotting the spike trains.
    h = findobj('tag', 'makeruns');
    if(isempty(h))
        h=figure('tag', 'makeruns');
    else
        figure(h);
        clf;
    end;
    for i = 1:ntestrep
        ya=i*ones(length(a.latency{i}));
        plot(a.latency{i}, ya, 'rd');
        hold on;
         yb=i*ones(length(b.latency{i}));
        plot(b.latency{i}, yb, 'b+');
    end;
    set(gca, 'Xlim', [0 1000]);
    set(gca, 'Ylim', [0 ntestrep+1]);
end;



% RANDPOIS(p,n)
%   This function creates a random variable with a Poisson distribution
%   with parameter 'a'.  If a second argument is used, a vector of
%   'n' poisson random variables is created.
%
% See also RAND, RANDN, RANDUNIFC, RANDEXPO, RANDGEO, RANDGAUSS

function out = randpois(a,n)

if nargin == 1
    i = 0;
    f = exp(-a);
    p = exp(-a);
    u = rand;
    while f<=u
        p = p * (a / (i + 1));
        f = f + p;
        i = i + 1;
    end
    out = i;
end

if nargin == 2
    for k=1:n
        i = 0;
        f = exp(-a);
        p = exp(-a);
        u = rand;
        while f<=u
            p = p * (a / (i + 1));
            f = f + p;
            i = i + 1;
        end
        out(k) = i;
    end
end

% RANDEXPO(a,n)
%   This function creates an exponentially distributed random variable
%   with parameter 'a'.  If a second argument is used, a vector of
%   'n' exponential variables is created.
%
% See also RAND, RANDN, RANDUNIFC, RANDGEO, RANDPOIS, RANDGAUSS

function out = randexpo(a, n)

if nargin == 1
    out = -log(rand) / a;
end

if nargin == 2
    randvec = rand(1,n);
    out = -log(randvec) ./ a;
end